<!DOCTYPE html>
<html lang="en">
  <head>
    <meta charset="utf-8">
    <title>Cindy JS</title>
    <script type="text/javascript" src="../../build/js/Cindy.js"></script>
    <script type="text/javascript" src="../../build/js/CindyGL.js"></script>
    <script type="text/javascript" src="../../build/js/symbolic.js"></script>
    <link rel="stylesheet" href="../../css/cindy.css">
  </head>
    
  <body style="font-family:Arial;">
  <h1>Raycasting Implicit Surfaces</h1>
  <div style="float: left; padding:20px">
      
      <script id="csmousedown" type="text/x-cindyscript">
      sx = mouse().x;
      sy = mouse().y;
      dragging = sx < .5;
      </script>
      <script id="csmouseup" type="text/x-cindyscript">
      dragging = false;
      </script>
      
      <script id="csinit" type="text/x-cindyscript">
      speciala = .5;
      smootha(x) := if(x<speciala,
        (-1/(speciala+.01))*(x-speciala)^2+speciala,
        (1/(1.01-speciala))*(x-speciala)^2+speciala      
      );
      
      seta(na) := (
        a = na;
        speciala = na;
        PA.y = (na-.5)*.7;
      );
      
      setzoom(zoom) := (
        PC.y = (zoom);
      );
      
      //initialize some variables
      mat = [
          [0.3513, -0.4908, -0.7973],
          [-0.8171, -0.5765, -0.0051],
          [-0.4571, 0.6533, -0.6036]
      ];
      sx = mouse().x;
      sy = mouse().y;
      dragging = false;
      N = 5;
      zoom = 2.2;
      a = 0.3;
      alpha = .7;
      
      //we stand at position mat*(0, 0, -2.2) and watch to (0,0,0).
      //ray(pixel, t) is the point in R^3 that lies at position t the ray behind the pixel at location pixel(vec2)
      //t=0 is corresponds to points within the interesting area near (0,0,0)
      ray(pixel, t) := mat * ((t+2.2) * (pixel.x, pixel.y, 1) + (0, 0, -2.2));
      
      //sphere with radius 1 for clipping
      S(r) := (r * r - 1);
      
      //fun is the user defined trivariate polynomial
      fun(x, y, z) := (x ^ 2 + y ^ 2 + z ^ 2 - (0.5 + a) ^ 2) ^ 2 - (3.0 * ((0.5 + a) ^ 2) - 1.0) / (3.0 - ((0.5 + a) ^ 2)) * (1 - z - sqrt(2) * x) * (1 - z + sqrt(2) * x) * (1 + z + sqrt(2) * y) * (1 + z - sqrt(2) * y);
      
      //F takes vec3 instead of 3 variables
      F(p) := (p=p/zoom;fun(p.x, p.y, p.z));
          
      //casteljau algorithm to evaluate and subdivide polynomials in Bernstein form.
      //poly is a vector containing the coefficients, i.e. p(x) = sum(0..N, i, poly_(i+1) * b_(i,N)(x)) where b_(i,N)(x) = choose(N, i)*x^i*(1-x)^(N-1)
      casteljau(poly, x) := (
        regional(alpha, beta);
        alpha = 1-x;
        beta = x;
        forall(0..N, k,
          repeat(N-k,
            poly_# = alpha*poly_# + beta*poly_(#+1);
          );
        );
        poly //the bernstein-coefficients of the polynomial in the interval [x,1]
      );

      //evaluates a polynomial, represented as vector of coefficients in bernstein-form
      eval(poly, x) := casteljau(poly, x)_1;
      
      //this function has to be called whenever fun changes
      init() := (
        dx = .05; dy =.02;
        diff(fun(x,y,z), x, dxF(x,y,z) := #);
        diff(fun(x,y,z), y, dyF(x,y,z) := #);
        diff(fun(x,y,z), z, dzF(x,y,z) := #);
        
          newN = degree(fun(x,y,z), x, y, z);
          if(newN==-1, newN=1000);
          if(newN!=N,
            N = newN;
            //The following line is DIRTY, but it makes the application run smooth for high degrees. :-)
            //Nethertheless, it might cause render errors for high degree surfaces. In fact, only a subset of the surface is rendered.
            //Adapt limit according to hardware.
            //values of kind 4*n-1 are good values, as it means to use vectors of length 4*n.
            N = min(N,7); 
            
            //N+1 Chebyshev nodes for interval (0, 1)
            li = apply(1..(N+1), k, (cos((2 * k - 1) / (2 * (N+1)) * pi)+1)/2);
            
            //A is the matrix of the linear map that evaluates a polynomial in bernstein-form at the Chebyshev nodes
            A = apply(li, node,
              //the i-th column contains the values of the (i,N) bernstein polynomial evaluated at the Chebyshev nodes
              apply(0..N, i, eval(
                apply(0..N, if(#==i,1,0)), // e_i = [0,0,0,1,0,0]
                node //evaluate  b_(i,N)(node)
              )) 
            );
            
            B = (inverse(A)); //B interpolates polynomials (in Bernstein basis), given the values [p(li_1), p(li_2), ...]
          )
          
      );
      init();
      
      //B3 is a matrix that interpolates quadratic polynomials (in monomial basis), given the values [p(-2), p(0), p(2)]
      B3 = inverse(apply([-2, 0, 2], c, apply(0 .. 2, i, c ^ i))); 

      //use symbolic differentation function
      dF(p) := (p=p/zoom; (
          dxF(p.x,p.y,p.z),
          dyF(p.x,p.y,p.z),
          dzF(p.x,p.y,p.z)
      ));

      //update the color color for the pixel at position pixel assuming that the surface has been intersected at ray(pixel, dst)
      //because of the alpha-transparency updatecolor should be called for the intersections with large dst first
      updatecolor(pixel, dst, color) := (
        regional(x, normal);
        x = ray(pixel, dst); //the intersection point in R^3
        color = (1 - alpha) * color;
              
        normal = dF(x);
        normal = normal / |normal|;

        forall(1..length(lightdirs),
          //illuminate if the normal and lightdir point in the same direction
          illumination = max(0, (lightdirs_# / abs(lightdirs_#)) * normal);
          color = color + alpha * (illumination ^ gamma_#) * colors_#;
        );
        color
      );
      

      nsign(pixel, a, b) := ( //Descartes rule of sign for the interval (a,b)
        //obtain the coefficients in bernstein basis of F along the ray in interval (a,b) by interpolation within this interval
        poly = B * apply(li,
          F(ray(pixel, a+#*(b-a))) //evaluate F(ray(pixel, ·)) along Chebyshev nodes for (a,b)
        );
        //count the number of sign changes
        ans = 0;
        //last = poly_1;
        forall(2..(N+1), k,
          //if(last == 0, last = poly_k;); this (almost) never happens
          if(min(poly_(k-1), poly_k) <= 0 & 0 <= max(poly_(k-1), poly_k), //sign switch; avoid products due numerics
            ans = ans + 1;
          );
        );
        ans //return value   
      );
      
      
      //bisect F(ray(pixel, ·)) in [x0, x1] assuming that F(ray(pixel, x0)) and F(ray(pixel, x1)) have opposite signs
      bisectf(pixel, x0, x1) := (
          regional(v0, v1, m, vm);
          v0 = F(ray(pixel, x0));
          v1 = F(ray(pixel, x1));
          repeat(11,
              m = (x0 + x1) / 2; vm = F(ray(pixel, m));
              if (min(v0,vm) <= 0 & 0 <= max(v0, vm), //sgn(v0)!=sgn(vm); avoid products due numerics
                  (x1 = m; v1 = vm;),
                  (x0 = m; v0 = vm;)
              );
          );
          m //return value   
      );
      
      //id encodes a node in a binary tree using heap-indices
      //1 is root node and node v has children 2*v and 2*v+1
      //computes s=2^depth of a node id: Compute floor(log_2(id));
      //purpose: id corresponds interval [id-s,id+1-s]/s
      gets(id) := (
        s = 1;
        repeat(15,
          if(2*s<=id,
            s = 2*s;
          )
        );
        s //return value
      );
      
      //determines the next node in the binary tree that would be visited by a regular in DFS
      //if the children of id are not supposed to be visited
      //In interval logic: finds the biggest unvisited interval directly right of the interval of id.
      next(id) := (
        id = id+1;
        //now: remove zeros from right (in binary representation) while(id&1) id=id>>1;
        repeat(15,
          if(mod(id,2)==0, 
            id = floor(id/2);
          )
        );
        if(id==1, 0, id) //return value - id 0 means we stop our DFS
      );
      </script>
      
      <script id="csdraw" type="text/x-cindyscript">
      //the following is executed for every rendered frame
      if (dragging,
          dx = 3 * (sx - mouse().x); dy = 3 * (sy - mouse().y);,
          dx = .9*dx; dy = .9*dy;
      );
      
      sx = mouse().x;
      sy = mouse().y;
      
      //the rotation matrix: It is modified either if the user is dragging or time passes
      mat = mat * (
          (1, 0, 0),
          (0, cos(dy), -sin(dy)),
          (0, sin(dy), cos(dy))
      ) * (
          (cos(dx), 0, -sin(dx)),
          (0, 1, 0),
          (sin(dx), 0, cos(dx))
      );
      
      
      //the 3 sliders at the left.
      PA.x = 0.55;
      if (PA.y > .4, PA.y = .4);
      if (PA.y < -.4, PA.y = -.4);
      a = smootha((.5 + PA.y/.7));

      PB.x = 0.6;
      if (PB.y > .4, PB.y = .4);
      if (PB.y < -.4, PB.y = -.4);
      alpha = ((.4 + PB.y) / .8) * .7 + .3;

      PC.x = 0.65;
      if (PC.y > .4, PC.y = .4);
      if (PC.y < -.4, PC.y = -.4);
      zoom = exp(3 * PC.y - 1);
      
      //configuration for the lights in the scene. A light has a position, a gamma-parameter for its shininess and a color
      lightdirs = [
          ray((.0, .0), -100), //enlights parts of the surface which normal points away from the camera
          ray((.0, .0), -100),
          ray((.0, .0), 100), //Has an effect, if the normal of the surface points to the camera
          ray((.0, .0), 100),
          (-10, 10, -2.),
          (10, -8, 3.)
      ];

      gamma = [1, 10, 1, 10, 5, 5];
      
              
      colors = [
          (.3, .5, 1.),
          (1, 2, 2) / 2,
          (1., 0.2, 0.1),
          (2, 2, 1) / 2,
          .9 * (.7, .8, .3),
          .9 * (.6, .1, .6)
      ];


      //what color should be given to pixel with pixel-coordinate pixel (vec2)
      if(alpha<.99,
        computeColor(pixel, l, u, color) := (
          regional(a, b);
          //traverse binary tree (DFS) using heap-indices
          //1 is root node and node v has children 2*v and 2*v+1
          id = 1; 
          //maximum number of steps
          repeat(min(newN*8,80),
            //id=0 means we are done; do only a DFS-step if we are not finished yet
            if(id>0,
              s = gets(id); //s = floor(log_2(id))
              
              //the intervals [a,b] are chossen such that (id in binary notation)
              //id = 1   => [a,b]=[l,u]
              //id = 10  => [a,b]=[l,(u+l)/2]
              //id = 101 => [a,b]=[l,(u+3*l)/4]
              //id = 11  => [a,b]=[(u+l)/2,u]
              //...
              a = u - (u-l)*((id+1)/s-1);
              b = u - (u-l)*((id+0)/s-1);
              
              //how many sign changes has F(ray(pixel, ·)) in (a,b)?
              cnt = nsign(pixel, a, b);
              if(cnt == 1 % (b-a)<1e-4, //in this case we found a root (or it is likely to have a multiple root)
                //=>colorize and break DFS
                color = updatecolor(pixel, bisectf(pixel, a, b), color);
                id = next(id),
              if(cnt == 0, //there is no root
                id = next(id), //break DFS
                
                //otherwise cnt>=2: there are cnt - 2*k roots.
                id = 2*id;  //visit first child within DFS
              )
          );  
          ));
          color
        );,
        computeColor(pixel, l, u, color) := (
          regional(a, b);
          //traverse binary tree (DFS) using heap-indices
          //1 is root node and node v has children 2*v and 2*v+1
          id = 1; 
          //maximum number of steps
          intersect = false;
          repeat(min(newN*7,50),
            //id=0 means we are done; do only a DFS-step if we are not finished yet
            if(!intersect & id>0,
              s = gets(id); //s = floor(log_2(id))
              
              //the intervals [a,b] are chossen such that (id in binary notation)
              //id = 1   => [a,b]=[l,u]
              //id = 10  => [a,b]=[l,(u+l)/2]
              //id = 101 => [a,b]=[l,(u+3*l)/4]
              //id = 11  => [a,b]=[(u+l)/2,u]
              //...
              a = l + (u-l)*((id+0)/s-1);
              b = l + (u-l)*((id+1)/s-1);
              
              //how many sign changes has F(ray(pixel, ·)) in (a,b)?
              cnt = nsign(pixel, a, b);
              if(cnt == 1 % (b-a)<1e-4, //in this case we found a root (or it is likely to have a multiple root)
                //=>colorize and break DFS
                intersect = true;
                id = next(id),
              if(cnt == 0, //there is no root
                id = next(id), //break DFS
                
                //otherwise cnt>=2: there are cnt - 2*k roots.
                id = 2*id;  //visit first child within DFS
              )
          );  
          ));
          if(intersect, updatecolor(pixel, bisectf(pixel, a, b), color), color)
        );
      );
        

      
      colorplot(
        spolyvalues = apply([-2, 0, 2], v, S(ray(#, v))); //evaluate S along ray
        spoly = B3 * spolyvalues;                         //interpolate to monomial basis
        D = (spoly_2 * spoly_2) - 4. * spoly_3 * spoly_1; //discriminant of spoly
        
        color = gray(0.7); //the color, which will be returned
        if (D >= 0, //ray intersects ball
          color = computeColor(
            #, 
            (-spoly_2 - re(sqrt(D))) / (2 * spoly_3), //intersection entering the ball
            (-spoly_2 + re(sqrt(D))) / (2 * spoly_3), //intersection leaving the ball
            color
          );              
        );
        color //return value
      ); //render the scene. # is the pixel coordinate
      
      
      drawtext((-.65, -.45), "degree: $" + if(newN<100,newN,"\infty") +"$");

      //lines for the sliders
      draw((.55, .4), (.55, -.4), color -> (0, 0, 0));
      draw((.6, .4), (.6, -.4), color -> (0, 0, 0));
      draw((.65, .4), (.65, -.4), color -> (0, 0, 0));    
      </script>

      <div  id="CSCanvas" style="border:0px"></div>
      
      <script type="text/javascript">
          var cdy = CindyJS({canvasname:"CSCanvas",
                      scripts: "cs*",
                      animation: {autoplay: true},
                      use : ["CindyGL","katex", "symbolic"],
                      geometry: [ { name:"PA", kind:"P", type:"Free", pos: [.5,.37,1], narrow: true, color: [1,1,1], size:8 },
                                  { name:"PB", kind:"P", type:"Free", pos: [.5,.5,1], narrow: true, color: [1,1,1], size:8 },
                                  { name:"PC", kind:"P", type:"Free", pos: [.5,.1,1], narrow: true, color: [1,1,1], size:8 } ],
                      ports: [{
                        id: "CSCanvas",
                        width: 700,
                        height: 500,
                        transform: [ { visibleRect: [ -0.7, -0.5, 0.7, 0.5 ] } ]
                      }]
          });
      </script>
      <input type="text" id="inp" value="(x^2+y^2+z^2-(0.5+a)^2)^2-(3*((0.5+a)^2)-1)/(3-((0.5+a)^2))*(1-z-sqrt(2)*x)*(1-z+sqrt(2)*x)*(1+z+sqrt(2)*y)*(1+z-sqrt(2)*y)"  onkeypress="if((event.which ? event.which : event.keyCode)==13) { cdy.evokeCS('fun(x,y,z) := (' + this.value + '); init();'); }" size="80" style="font-size:18px">
      </div>
  <p>Non-algebraic functions are approximated by polynomials.
  Roots are isolated by Descartes Method in Bernstein basis.</p>
  <p>You can enter your own implicit surfaces or select one of the following list:</p>
  <p><select id="sel" size="15" style="width:20em;"><option data-a="1" value="(x^2+y^2+z^2-(0.5+a)^2)^2-(3*((0.5+a)^2)-1)/(3-((0.5+a)^2))*(1-z-sqrt(2)*x)*(1-z+sqrt(2)*x)*(1+z+sqrt(2)*y)*(1+z-sqrt(2)*y)">Kummer Quartic</option>
    <option data-a="1" value="4*((a*(1+sqrt(5))/2)^2*x^2-y^2)*((a*(1+sqrt(5))/2)^2*y^2-z^2)*((a*(1+sqrt(5))/2)^2*z^2-x^2)-1*(1+2*(a*(1+sqrt(5))/2))*(x^2+y^2+z^2-1)^2">Barth Sextic</option>
    <option data-a="0" value="-2*a/125+x^8+y^8+z^8-2*x^6-2*y^6-2*z^6+1.25*x^4+1.25*y^4+1.25*z^4-0.25*x^2-0.25*y^2-0.25*z^2+0.03125">Chmutov Octic</option>
    <option data-a="1" data-zoom="-.1" value="a*(-1/4*(1-sqrt(2))*(x^2+y^2)^2+(x^2+y^2)*((1-1/sqrt(2))*z^2+1/8*(2-7*sqrt(2)))-z^4+(0.5+sqrt(2))*z^2-1/16*(1-12*sqrt(2)))^2-(cos(0*pi/4)*x+sin(0*pi/4)*y-1)*(cos(pi/4)*x+sin(pi/4)*y-1)*(cos(2*pi/4)*x+sin(2*pi/4)*y-1)*(cos(3*pi/4)*x+sin(3*pi/4)*y-1)*(cos(4*pi/4)*x+sin(4*pi/4)*y-1)*(cos(5*pi/4)*x+sin(5*pi/4)*y-1)*(cos(6*pi/4)*x+sin(6*pi/4)*y-1)*(cos(7*pi/4)*x+sin(7*pi/4)*y-1)">
      Endraß Octic
    </option>
    <option data-zoom=".2" value="x^2+y^2+z^2-1">Ball</option>
    <option data-zoom=".2" value="k = 6; x^k+y^k+z^k-1">Cube</option>
    <option data-zoom=".2" value="x^2+z^2-1/3*(1+y)^3*(1-y)^3">Citric</option>
    <option data-zoom=".1" value="x^2+y^2+z^3-z^2">Ding Dong</option>
    <option data-zoom="0" value="x^3+x^2*z^2-y^2">Hummingbird</option>
    <option data-zoom=".2" value="x^2-x^3+y^2+y^4+z^3-z^4">Vis a Vis</option>
    <option data-zoom=".1" value="(x^2+9/4*y^2+z^2-1)^3-x^2*z^3-9/80*y^2*z^3">Sweet</option>
    <option data-zoom=".2" data-a="1/4" value="k=a*2;(x+(k/2-1))*(x^2+y^2+z^2-k^2/4)+z^2">Parabolic Ring Cyclide</option>
    <option data-a="0" data-zoom="-.15" value="cos(x)*sin(y) + cos(y)*sin(z) + cos(z)*sin(x) + a">Gyroid</option>
    <option data-a="0" data-zoom="-.15" value="cos(x)+cos(y)+cos(z)+a">Schwarz P</option>
    <option data-a=".5" data-zoom=".1"  value="r=a; R=1; (x^2+y^2+z^2+R^2-r^2)^2-4*R^2*(x^2+y^2)">Torus</option>
    <option data-a=".4" data-zoom="-.1" value = "r=a/2; R=.9; ((sin(x)^2+y^2+z^2+R^2-r^2)^2-4*R^2*(sin(x)^2+y^2))*((cos(x)^2+y^2+z^2+R^2-r^2)^2-4*R^2*(cos(x)^2+z^2))">Interleaved Tori</option>
  </select></p>
  <script type="text/javascript">
  var select = document.getElementById("sel");
  select.addEventListener('change', function(event) {
    document.getElementById('inp').value = this.value;
    cdy.evokeCS('fun(x,y,z) := (' + this.value + '); init();');

    var a = this.options[this.selectedIndex].getAttribute("data-a") || .5;
    cdy.evokeCS('seta(' + a + ')');

    var zoom = this.options[this.selectedIndex].getAttribute("data-zoom");
    if(zoom) cdy.evokeCS('setzoom(' + zoom + ')');

  }, false);
  </script>

  </body>
</html>
